Computation of a STAI dataset with Field II to simulate specular tissue using parameters of an L11-4v 128 element Verasonics Transducer and beamforming with USTB and using MATLAB's parfor function for parallelization
This example shows how to load the data from a Field II simulation of specular tissue into USTB objects, and then beamform it with the USTB routines. This example uses the L11-4v 128 element Verasonics Transducer The Field II simulation program (field-ii.dk) should be in MATLAB's path.
This tutorial assumes familiarity with the contents of the 'CPWC simulation with the USTB built-in Fresnel simulator' tutorial. Please feel free to refer back to that for more details.
by Alfonso Rodriguez-Molares alfonso.r.molares@ntnu.no, Ole Marius Hoel Rindal olemarius@olemarius.net and Arun Asokan Nair anair8@jhu.edu 09.05.2017
Contents
- Close old plots
- Basic Constants
- field II initialisation
- Transducer definition L11-4v, 128-element linear array transducer
- Pulse definition
- Aperture Objects
- Speckle Phantom
- Output data
- Compute STA signals
- Channel Data
- Save UFF dataset
- Scan
- Beamformer
- Plot a PDF of speckle and compare to theoretically predicted Rayleigh distribution that speckle possesses
Close old plots
close all; filename = 'FieldII_speckle_simulation.uff'; if exist([data_path filesep filename]) > 0 %If the file exists load the file disp('We were lucky, and the data was already simulated, so we can simply load it!'); uff_file= [data_path filesep filename]; channel_data = uff.read_object(uff_file,'/channel_data'); else % Else, run the simulation
Basic Constants
Our first step is to define some basic constants for our imaging scenario - below, we set the speed of sound in the tissue, sampling frequency and sampling step size in time.
c0=1540; % Speed of sound [m/s] fs=100e6; % Sampling frequency [Hz] dt=1/fs; % Sampling step [s]
field II initialisation
Next, we initialize the field II toolbox. Again, this only works if the Field II simulation program (field-ii.dk) is in MATLAB's path. We also pass our set constants to it.
field_init(0);
set_field('c',c0); % Speed of sound [m/s]
set_field('fs',fs); % Sampling frequency [Hz]
set_field('use_rectangles',1); % use rectangular elements
*------------------------------------------------------------*
* *
* F I E L D I I *
* *
* Simulator for ultrasound systems *
* *
* Copyright by Joergen Arendt Jensen *
* Version 3.30, April 5, 2021 (Matlab 2021a version) *
* Web-site: field-ii.dk *
* *
* This is citationware. Note the terms and conditions *
* for use on the web-site at: *
* field-ii.dk/?copyright.html *
* It is illegal to use this program, if the rules in the *
* copyright statement is not followed. *
*------------------------------------------------------------*
Warning: Remember to set all pulses in apertures for the new sampling frequency
Transducer definition L11-4v, 128-element linear array transducer
Our next step is to define the ultrasound transducer array we are using. For this experiment, we shall use the L11-4v 128 element Verasonics Transducer and set our parameters to match it.
probe = uff.linear_array();
f0 = 5.1333e+06; % Transducer center frequency [Hz]
lambda = c0/f0; % Wavelength [m]
probe.element_height = 5e-3; % Height of element [m]
probe.pitch = 0.300e-3; % probe.pitch [m]
kerf = 0.03e-03; % gap between elements [m]
probe.element_width = probe.pitch-kerf;% Width of element [m]
lens_el = 20e-3; % position of the elevation focus
probe.N = 128; % Number of elements
pulse_duration = 2.5; % pulse duration [cycles]
Pulse definition
We then define the pulse-echo signal which is done here using the fresnel simulator's pulse structure. We could also use 'Field II' for a more accurate model.
pulse = uff.pulse();
pulse.center_frequency = f0;
pulse.fractional_bandwidth = 0.65; % probe bandwidth [1]
t0 = (-1/pulse.fractional_bandwidth/f0): dt : (1/pulse.fractional_bandwidth/f0);
impulse_response = gauspuls(t0, f0, pulse.fractional_bandwidth);
impulse_response = impulse_response-mean(impulse_response); % To get rid of DC
te = (-pulse_duration/2/f0): dt : (pulse_duration/2/f0);
excitation = square(2*pi*f0*te+pi/2);
one_way_ir = conv(impulse_response,excitation);
two_way_ir = conv(one_way_ir,impulse_response);
lag = length(two_way_ir)/2;
% We display the pulse to check that the lag estimation is on place
% (and that the pulse is symmetric)
figure;
plot((0:(length(two_way_ir)-1))*dt -lag*dt,two_way_ir); hold on; grid on; axis tight
plot((0:(length(two_way_ir)-1))*dt -lag*dt,abs(hilbert(two_way_ir)),'r')
plot([0 0],[min(two_way_ir) max(two_way_ir)],'g');
legend('2-ways pulse','Envelope','Estimated lag');
title('2-ways impulse response Field II');
Aperture Objects
Next, we define the the mesh geometry with the help of Field II's xdc_linear_array function.
noSubAz=round(probe.element_width/(lambda/8)); % number of subelements in the azimuth direction noSubEl=round(probe.element_height/(lambda/8)); % number of subelements in the elevation direction Th = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]); Rh = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]); % We also set the excitation, impulse response and baffle as below: xdc_excitation (Th, excitation); xdc_impulse (Th, impulse_response); xdc_baffle(Th, 0); xdc_center_focus(Th,[0 0 0]); xdc_impulse (Rh, impulse_response); xdc_baffle(Rh, 0); xdc_center_focus(Rh,[0 0 0]);
Speckle Phantom
In our next step, we define our phantom. Here, our goal is to simulate speckle so we have a 100 scatterrers with axial and lateral coordinates randomly drawn from a uniform distribution and scatterer amplitudes randomly drawn from a normal distribution.
number_of_scatterers = 500000;
xxp_speckle=random('unif',-5e-3,5e-3,number_of_scatterers,1);
zzp_speckle=random('unif',15e-3,20e-3,number_of_scatterers,1);
sca = [xxp_speckle zeros(length(xxp_speckle),1) zzp_speckle]; % list with the scatterers coordinates [m]
amp=randn(length(sca),1); % list with the scatterers amplitudes
cropat=round(1.1*2*sqrt((max(sca(:,1))-min(probe.x))^2+max(sca(:,3))^2)/c0/dt); % maximum time sample, samples after this will be dumped
Output data
We define the variables to store our output data
t_out=0:dt:((cropat-1)*dt); % output time vector STA=zeros(cropat,probe.N,probe.N); % impulse response channel data
Compute STA signals
Now, we finally reach the stage where we generate a STA (Synthetic Transmit Aperture) dataset with the help of Field II.
disp('Field II: Computing STA dataset'); disp('No waitbar possible for parfor, so just be patient :)'); parfor n=1:probe.N %Since we are using parfor, we have to initate Field II and the arrays %for every worker as well. field_init(0); Th = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]); Rh = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]); xdc_excitation (Th, excitation); xdc_impulse (Th, impulse_response); xdc_baffle(Th, 0); xdc_center_focus(Th,[0 0 0]); xdc_impulse (Rh, impulse_response); xdc_baffle(Rh, 0); xdc_center_focus(Rh,[0 0 0]); % transmit aperture xdc_apodization(Th, 0, [zeros(1,n-1) 1 zeros(1,probe.N-n)]); xdc_focus_times(Th, 0, zeros(1,probe.N)); % receive aperture xdc_apodization(Rh, 0, ones(1,probe.N)); xdc_focus_times(Rh, 0, zeros(1,probe.N)); % do calculation [v,t]=calc_scat_multi(Th, Rh, sca, amp); % save data -> with parloop we need to pad the data if size(v,1)<cropat STA(:,:,n)=padarray(v,[cropat-size(v,1) 0],0,'post'); else STA(:,:,n)=v(1:cropat,:); end % Sequence generation seq(n)=uff.wave(); seq(n).probe=probe; seq(n).source.xyz=[probe.x(n) probe.y(n) probe.z(n)]; seq(n).sound_speed=c0; seq(n).delay = probe.r(n)/c0-lag*dt+t; % t0 and center of pulse compensation end
Index exceeds the number of array elements (1).
Error in STAI_L11_speckle_parfor (line 139)
disp('Field II: Computing STA dataset');
Channel Data
In this part of the code, we creat a uff data structure to specifically store the captured ultrasound channel data.
channel_data = uff.channel_data();
channel_data.sampling_frequency = fs;
channel_data.sound_speed = c0;
channel_data.initial_time = 0;
channel_data.pulse = pulse;
channel_data.probe = probe;
channel_data.sequence = seq;
channel_data.data = STA./max(STA(:));
Save UFF dataset
Finally, we save the data into a UFF file. There is a
channel_data.write([data_path filesep filename],'channel_data');
end
Scan
The scan area is defines as a collection of pixels spanning our region of interest. For our example here, we use the linear_scan structure, which is defined with two components: the lateral range and the depth range. scan too has a useful plot method it can call.
scan=uff.linear_scan('x_axis',linspace(-5e-3,5e-3,256).', 'z_axis', linspace(15e-3,20e-3,256).');
Beamformer
With channel_data and a scan we have all we need to produce an ultrasound image. We now use a USTB structure beamformer, that takes an apodization structure in addition to the channel_data and scan.
pipe=pipeline(); pipe.channel_data=channel_data; pipe.scan=scan; % Delay and sum on receive, then coherent compounding b_data=pipe.go({midprocess.das() postprocess.coherent_compounding()}); % Display image figure(1);clf b_data.plot(1)
Plot a PDF of speckle and compare to theoretically predicted Rayleigh distribution that speckle possesses
envelope = abs(b_data.data); envelope = envelope./max(envelope(:)); m = mean(envelope(:)); s = std(envelope(:)); snr_calculated_das = m/s snr_theoretical = (pi/(4-pi))^(1/2) b = s/(sqrt((4-pi)/2)); %Scale parameter % Estimate PDF x_axis = linspace(0,1,200); [n,xout] = hist(envelope(:),x_axis); delta_x = xout(2)-xout(1); n = n/sum(n)/delta_x; % Theoretical Rayleigh PDF theoretical_pdf = (x_axis./b^2).*exp(-x_axis.^2/(2.*b^2)); % Plot color=[0.25 1 0.75] figure(2);clf; plot(xout,n,'LineWidth',2,'Color','r','DisplayName','Estimated PDF');hold on; plot(x_axis,theoretical_pdf,'--','Color',color,'LineWidth',2,'DisplayName','Rayleigh Theoretical PDF'); title('PDF of envelope'); xlabel('Normalized amplitude'); ylabel('Probability') legend('show');